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Abstract 

Molecular dynamics simulation is used to study the time-scales involved in the homoge- 
neous melting of a superheated crystal. The interaction model used is an embedded-atom 
model for Fe developed in previous work, and the melting process is simulated in the 
microcanonical (N, V, E) ensemble. We study periodically repeated systems containing 
from 96 to 7776 atoms, and the initial system is always the perfect crystal without free 
surfaces or other defects. For each chosen total energy E and number of atoms N, we per- 
form several hundred statistically independent simulations, with each simulation lasting 
for between 500 ps and 10 ns, in order to gather statistics for the waiting time r w before 
melting occurs. We find that the probability distribution of t w is roughly exponential, and 
that the mean value (r w ) depends strongly on the excess of the initial steady temperature 
of the crystal above the superheating limit identified by other researchers. The mean (r w ) 
also depends strongly on system size in a way that we have quantified. For very small 
systems of ~ 100 atoms, we observe a persistent alternation between the solid and liquid 
states, and we explain why this happens. Our results allow us to draw conclusions about 
the reliability of the recently proposed Z method for determining the melting properties 
of simulated materials, and to suggest ways of correcting for the errors of the method. 

1 Introduction 

The supercooling of liquids below their thermodynamic freezing point is familiar and easily 
observable, but the superheating of solids above their melting point is much more difficult 
to study. This is because melting is usually initiated at surfaces (grain boundaries and 
other defects may also initiate melting) , so that superheating is generally possible only in 
solids that have no surfaces [T] . Melting from the defect-free superheated state, sometimes 
called "homogeneous melting", has been experimentally observed 2, 3, HIE], but there is 
still rather little detailed understanding of the kinetics of the process. Fortunately, com- 
puter simulation offers a rather straightforward way of studying superheated defect-free 
solids, and this has led to a recent resurgence of interest in the subject [B]. In addition to 
the purely scientific interest, it has recently been shown that the concept of the "super- 
heating limit" leads naturally to a simulation technique known as the "Z method" that 
offers a new and potentially useful way of determining the melting properties of simulated 



materials [7]. In the present paper, we report new simulation results on the kinetics of 
homogeneous melting which shed light on the conditions needed for the Z method to yield 
reliable results. 

The study of melting properties by computer simulation dates back over 50 years [5J. 
Two main approaches have become firmly established over that period. The first relies on 
the separate calculation of the free energies of the solid and liquid, and is based on the fact 
that the two phases coexist in thermal equilibrium when both the pressures and the chem- 
ical potentials in the solid and liquid are equal [HI EH E] ■ The second approach consists 
of the explicit simulation of coexisting solid and liquid in the same simulated system j!2j. 
For any given interaction model, careful application of the two approaches to large enough 
systems should yield essentially identical results for the relation between melting tempera- 
ture T m and the pressure P on the melting curve, as well as other properties, including the 
heat and volume of fusion. The two approaches have been extensively used to determine 
the melting properties of a wide variety of systems interacting via classical potentials, 
including hard spheres [5J [T3] , inverse-power [TJJ [TS] and Lennard- Jones models [TB] , as 
well as more complex models such as the Born-Mayer model of ionic liquids |17j , a variety 
of models for water [18], and the embedded-atom model for metals 19, 20 . In the past 
15 years, there has been rapidly increasing interest in the determination of melting prop- 
erties using ab initio molecular dynamics simulation (AIMD) based on density-functional 
theory (DFT) . Initially, the free-energy route was used [21] [251 H3] > with thermodynamic 
integration employed to compute the difference of free energy between the ab initio system 
and an appropriately chosen reference system. The coexistence approach has also been 
extensively used, mainly with paramcteriscd classical potentials tuned to data produced 
by DFT simulations on the solid and liquid. However, there have also been several studies 
in which direct ab initio simulations of coexisting solid and liquid have been performed 
on systems of several hundred atoms [231 [231 1251 [271 I28"II29] . 

The point of departure of the Z method 7 is an apparently simple question about 
the superheating of a solid: If a solid is allowed to evolve under the classical equations 
of motion at constant number of atoms TV, volume V and internal energy E, what is 
the maximum energy Tvls it can have without eventually transforming completely into 
the liquid state? The proposed answer [7] is that it is the lowest energy on the given 
isochorc within the field of thermodynamic stability of the liquid [30 . This energy E^g 
corresponds to a temperature Tls representing the limit of superheating, above which the 
solid, evolving at constant energy, will always melt. Since £ls is the lowest energy of the 
liquid on the isochore, it should be the energy of the liquid in coexistence with the solid, 
so that the energy Tls of the liquid is associated with the melting temperature T m : 

E sol (V,T LS ) = E h *(V,T m ) . (1) 

The procedure for determing the point (T m , P) on the melting curve belonging to a speci- 
fied liquid-state V is then to perform a sequence of (N, V, E) m.d. simulations, monitoring 
T and P in each, the aim being to locate the threshold Tls (equivalently, the threshold 
Tls), above which the transition to the liquid always occurs, and below which it never 
occurs [7j. 

Implicit in these statements is an important question about timescales, i.e. about the 
kinetics of homogeneous melting. The original papers about the Z method [7] emphasise 
that in order to be reasonably certain whether the initial T is above or below Tls > evolution 
must be allowed over a long enough time, which may be on the order of ns, and the number 
of atoms N must also be large enough. (For simplicity, we assume here a single-component 
system of atoms.) This naturally raises a number of important questions, which we shall 
try to answer. First, since homogeneous melting in constant- (TV, V, E) dynamics appears 
to be a rare-event process, we want to examine the probability distribution of waiting 
times r w before the transition occurs. This means repeating the simulation many times 
with the same (JV, V, E) but with statistically independent initial conditions. Second, 
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we want to study how the mean waiting time (r w ) depends on how far above Tls the 
system is initiated. Third, we need to determine the dependence of (t w ) on the system 
size N. Naturally, the numerical answers to these questions will depend on the nature 
of the system and the number density n = N/V. Given the recent interest in using the 
Z method for the melting of metals [31] , particularly at high pressures [32J 133] > we have 
decided to study the statistics of r w using an embedded atom model (EAM) for Fe, whose 
melting properties are already well known from previous work |34j . We will also present 
Z-mcthod calculations on the melting of Fe using AIMD. We take a density corresponding 
to the megabar pressures typical of the Earth's core. 

In the next Section, we summarise the details of the interaction model and the simu- 
lation procedures. Our results on the statistical distribution of r w and the dependence of 
(r w ) on initial and final temperatures and system size from EAM and AIMD simulations 
are presented in Sec. 3. In the final Section, we discuss the implications of our results 
for the understanding of homogeneous melting and for the reliability of the Z method, 
particularly in the context of ab initio simulations. 

2 Techniques 

The embedded-atom model (EAM) for Fe used in our simulations is the one used as 
a reference system in our earlier work [34 on the ab initio melting curve of hep Fe. 
Essentially the same EAM was also used in the very recent work of Belonoshko et al. [33], 
in which they used the Z method to study the melting of Fe and an Fe/Si alloy. The model 
is actually a modification of a much earlier EAM developed originally by Belonoshko's 
group [35] . We recall that in an EAM scheme the total potential energy E tot is expressed 
as a sum of atomic energies: -Etot = X)j ^ii with the sum running over the N atoms in 
the system. Each term is a sum of two parts: Ei = E™ p + F{pi). Here, El cp consists of 
a sum of repulsive inverse-power pair potentials: E\ cp = X)j e ( a / r 'j.j) n > where is the 
distance between atoms % and j, and the term i = j is excluded. F(pi) is an "embedding 
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function" which describes the metallic bonding. It has the form F(pi) — —eCp t , with 
Pi = J2j' (a/rij)" 1 . The values of a and m are those in the original Belonoshko model [35] . 
while in Ref. |34j we showed how all the other parameters could be optimised by minimising 
the fluctuations of the difference between the EAM and ab initio energies in simulations 
of the liquid and the high-temperature solid. The numerical values of the parameters are: 
a = 3.4714 A, m = 4.788, e = 0.1662 eV, n = 5.93, C = 16.55. We apply the spatial 
cut-off r c — 5.5 A, so that terms in both £ , l rop and pi for which r,j > r c are set to zero, 
with the usual cutting and shifting used to ensure continuity. 

Following previous work [35] [33 [3H] , we focus here on the melting properties of Fe in 
the high-pressure region that is important for understanding the solid inner core and the 
liquid outer core of the Earth. Specifically, we confine ourselves to pressures P ~ 330 GPa, 
which is the pressure at the boundary between the inner and outer core |39j . From 
extensive simulations with our EAM on large systems containing solid and liquid in stable 
coexistence [31], we know that its melting temperature at P = 323 GPa is 6200 ± 100 K. 
We have recently refined these coexistence simulations so as to reduce the statistical errors, 
finding that a more accurate value of T m at this pressure is 6215 K. 

All the simulations to be presented were performed at the same density corresponding 
to a volume per atom V/N = 7.139 A 3 , which gives pressures in the region of interest. In 
every simulation, we start from the perfect hep crystal, with all atoms on their regular 
lattice sites, and we assign random velocities drawn from a Maxwellian distribution, the 
velocities then being shifted and scaled so that the total momentum is zero and the kinetic 
energy per atom K/N has a value corresponding exactly to a chosen initial temperature 
2] = 2K/3k#N . Verlet's algorithm [JU] was used with a time-step of 1 fs, which ensures 
conservation of total energy with a drift of typically no more than ~ 10 K over times of 
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several ns. 

We shall present simulations on systems containing N = 96, 150, 392, 972 and 7776 
atoms. For each N and for each initial temperature T iy we have performed several hundred 
m.d. simulations of at least 500 ps (in some cases we have continued the simulations 
for over 10 ns), in order to gather statistical information about the melting process. 
For the larger systems, the simulations were run on large parallel computers, with each 
individual simulation running on 24 cores, and with typically 350 such simulations running 
simultaneously. This mode of operation makes it possible, for example, to run an overall 
total of ~ 1 [is of m.d. on the 7776-atom system in only a few hours of wall-clock time. 
For the small systems, we ran the calculations on single processors using local facilities. 

The ab-initio simulations were run with exactly the same technical details as described 
in [36l[37l|38], using the VASP code [41] with the projector augmented wave method [421143] 
and an efficient extrapolation of the charge density [44] . 

3 Results 

We start by showing some examples of homogeneous melting from our simulations on 
the system of 7776 atoms. Fig. [I] displays the time-dependent temperature and pressure 
in four simulations that were all initiated from the perfect hep crystal with exactly the 
same kinetic energy corresponding to the temperature 7] = 15600 K, but with statistically 
independent random velocities. As expected from equipartition, T drops rapidly to about 
half its initial value (this rapid drop is not shown in the Figure), and it then fluctuates 
about a quasi-steady value T so \ = 7590 K for many ps, before it drops again over a rather 
short period of ~ 8 ps, and then fluctuates again about a lower steady value Tu q = 6315 K. 
The second drop is due to melting, as we have verified by monitoring the self-diffusion 
of atoms via the time-dependent mean-square displacement (Ar(t) 2 ). The appearance of 
atomic disorder throughout the system when the system melts is also easy to observe in 
movies prepared from the coordinate files. Melting is accompanied by an increase of P 
by ~ 10 GPa, which occurs over the same rather short interval as the drop in T. These 
effects are familiar from many previous reports on the Z method [71 1321 133) : the drop in 
T is due to the latent heat of fusion, and the increase of P is associated with the volume 
increase that would occur on melting if the pressure were held constant. We note that the 
waiting times r w that elapse before melting are different in the four examples shown. This 
is what we expect of a rare-event process, and is consistent with the statements in earlier 
reports [7J [321 133] that the time at which the melting transition occurs is not correlated 
with the details of the initial conditions. We find that the final mean temperature Tn q 
and pressure Pu q of the liquid are the same in all the examples. This is expected, because 
in every case the system settles into exactly the same thermodynamic state of the liquid. 
The temperature Tu q is somewhat above the melting temperature T m at pressure P\i q , as 
expected because the system was initiated above the limit of superheating. 

These observations naturally raise the question of the statistical distribution of waiting 
times r w . To investigate this, we have to repeat the simulations many times, starting al- 
ways from the perfect lattice with exactly the same initial kinetic energy, but independent 
random velocities. To make the question well posed, we need a definition of r w . We note 
from Fig.[T]that the fluctuations of T about its mean value in the quasi-steady state of the 
solid before melting are much smaller than the drop of T during the melting process. We 
therefore define t w to be the elapsed time from the start of the simulation to the instant 
when T is mid-way between the mean quasi-steady temperature T so \ of the solid and the 
mean final temperature Ti; q of the liquid. 

For the 7776-atom system with 7] ~ 15800 K, we repeated the simulations 350 times, 
with each run having a duration of 650 ps. We found that it melted in all cases, and we 
accumulated the histogram of t w shown in Fig.[2j We see that after a short incubation time 
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of no more than ~ 20 ps, the probability distribution of r w decays in a quasi-exponential 
way. This is what we expect if melting is a random process having short memory time 
with a constant probability per unit time 1 /tq of occurring, given that it has not already 
occurred. In this case, the probability distribution of waiting times p(r w ) would have the 
form Tq 1 exp(— t w /to) and the mean waiting time would be (r w ) = To. We show in Fig[2] 
a fit of the exponential function to the histogram. The value of To given by this fit is 
To = 24.1 ps, which agrees well with the value (t w ) = 24.7 ps computed directly from 
the sample of 350 values of t w . We have checked that the final mean 7]i q and Pu q of 
the liquid are the same in all cases, within statistical error, having the numerical values 
Tn q = 6410 ± 5 K and Pi = 330.4 ± 0.3 GPa. 

For comparison, we show the results for p(r w ) when we do exactly the same thing 
for the 7776-atom system, but now with a higher initial T| of 16000 K, the number of 
independent simulations in this case also being 350. As expected, melting occurs on 
average more rapidly with this Tj, the values of To from the exponential fit and from the 
directly computed (t w ) being 8.3 ps and 9.4 ps. The mean temperatures of the solid and 
the liquid in this case are are T so i = 7750 K and 7ii q = 6510 K, which, as expected, are 
higher than the values found with the lower I). 

It is clear from these observations that the mean waiting time depends on the extent 
to which the temperature T so i exceeds the limit of superheating Tls. To study this 
further, we have repeated the simulations with 7] values of 16000, 15800, 15600, 15400 
and 15200 K. At the lowest of these 7], melting was not seen in any of the simulations, even 
though we repeated them 350 times with statistically independent initial velocities, the 
duration of the simulation being 710 ps in every case. At 7] = 15400 K and 7] = 15600 K , 
melting was observed within 660 ps in only 14 and 283 out of 350 simulations, respectively. 
At all the other 7] values, melting occurred in all cases, and we were able to construct 
essentially complete histograms; the values of (t w ) and the value of To obtained by fitting 
p(t w ) = T^ 1 exp(— </t ) to the histogram agreed closely. 

The final mean 7]i q of the liquid is a monotonically increasing function of 7], and 
it is convenient to examine the dependence of (t w ) on 71i q . We have found that it is 
helpful to plot the quantity (t w ) -1 / 2 against 7]i q , as shown in Fig. k5[ The points fall 
roughly on a straight line, and the indication is that (t w ) -1 / 2 — > (lTe. (t w ) — > oo) at 
a characteristic temperature. At the same time, Pn q also tends to a limiting value. We 
identify the characteristic temperature as the melting temperature T m at the pressure 
Py lq , because T m is the lowest possible final mean temperature of the liquid, namely the 
temperature found when T so i = T^g. The value of T m obtained by this extrapolation is 
6260 K, the extrapolated pressure being 328 GPa. These results agree very well with the 
value T m = 6215 K from explicit coexistence simulations at the pressure P = 323 GPa 
(see Sec. [2]). 

All the results presented so far are for the large system of 7776 atoms. We have 
performed simulations of essentially the same kind for systems of N = 972, 392, 150 and 
96 atoms, in each case initiating the simulations at sequence of initial temperatures 7], 
repeating the simulations at each (N, 7]) a few hundred times, determining the liquid-state 
7ii q values for the cases where melting occurs, and extracting the (rJ) values. The plots 
of (t w ) -1 / 2 against Tn q for all the system sizes are displayed in Fig. [3J The results appear 
to be very coherent: for each N, the (t w ) -1 / 2 points fall reasonably well on a straight line, 
and the Figure shows the linear least-square fits. These linear fits extrapolate to give T m 
values that agree for the different system sizes to within ~ 100 K, i.e. to within ~ 2 %. 
This agreement suggests that the Z method can be a robust way of obtaining T m values 
close to the thermodynamic limit for systems that would be much too small for explicit 
coexistence simulations, provided very long simulations are performed and provided one 
extrapolates to the (t w ) — > 00 limit. 

The plots of Fig. |3] show that for the given density n the mean waiting times (t w ) for 
different degrees of superheating and different system sizes can all be roughly represented 
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by the formula (r w ) = A/(Ty lq — T m ) 2 , where T m is independent of Ty lq and N but A 
depends on N. They also show that A increases with decreasing N. More extensive 
results would be needed to make precise statements about this N dependence, but we 
find that the inverse proportionality A = B/N fits the results quite well. As evidence for 
this, we display in Fig. [4 a plot of (iV^w))" 1 / 2 against Tu q , showing that all our data are 
quite well reproduced by the formula: 

(A^Tw))- 1 / 2 = C(T liq -T m ) , (2) 

where T m = 6315 K and C = B^ 1 / 2 has the value 1.9 x 1CT 5 ps" 1 / 2 K _1 . 

For the system of 96 atoms, we observe a significant and interesting effect, which 
sheds further light on the kinetics of homogeneous melting. After the system has made 
the transition from superheated solid to liquid, it remains in the liquid state only for a 
finite time, before reverting back to the solid state. In fact, if the simulation is continued 
long enough, we observe a continual alternation between the solid and liquid states. An 
example of this behaviour is shown in Fig. [5] where we see that the temperature alternates 
between solid- like and liquid-like values T so \ and Iii q . This effect becomes very clear if 
we construct a histogram of temperature T obtained by sampling over the course of 
many simulations, all having exactly the same total energy E, and hence the same liquid 
temperature Tn q . The temperature histograms for the 96-atom system for different Tn q 
values are shown in Fig.[6j We see that in each case T has a bimodal distribution, being the 
superposition of the quasi-Gaussian distributions that would be found if the system were 
wholly in the solid state or wholly in the liquid state. In fact, we can fit the histograms 
very well by a superposition of Gaussians, and the relative weights of the two Gaussians 
for a given Tn q tell us the relative amounts of time spent by the system in the solid and 
liquid states. We display the fraction of time ai; q spent in the liquid state as a function 
of T liq in Fig. [7) 

At first sight, the ceaseless alternation between solid and liquid might seem surpris- 
ing, because it implies that homogeneous melting from the superheated solid is not the 
irreversible process that one might expect. However, in Sec. [4] we will point out why this 
alternation is required by the principles of statistical mechanics, we give a simple formula 
that explains why ai; q depends on Xi iq as shown in Fig. [7J and we discuss whether the 
alternation should also be seen in larger systems. 

There has been considerable interest in using the Z-method with AIMD simulation 
(we refer to this as AIMD-Z) to obtain melting curves, particularly for metals [32 , 33J. To 
test the practical operation of AIMD-Z, we have performed our own calculations on the 
high-P melting of hep Fe, on which there is already very extensive previous work based 
on both free-energy and explicit-coexistence methods, including a recent study of AIMD 
coexistence on systems of N — 980 atoms [28 . Before presenting our AIMD-Z results 
on this problem, it is useful to consider the errors that can be expected. Our AIMD-Z 
simulations were performed on systems of 150 atoms with duration of ~ 50 ps. Clearly, 
with a run of this duration, initiated above T^g, we are unlikely to observe homogeneous 
melting unless (r w ) is ~ 50 ps or less. From the formula given in eqn ([2]), we estimate that 
this will yield Tn q — T m ~ 600 K, and this is the error we may make if T m is estimated as 
the lowest Ti; q for which homogeneous melting is observed. 

We present in Fig.|H]our AIMD-Z results for the melting of hep Fe with N = 150 using 
50-ps simulations. The Figure also shows the melting curve obtained many years ago with 
exactly the same AIMD techniques but based on free-energy calculations [37J 13H1 El] > as 
well as a point on the melting curve from AIMD coexistence using 980 atoms 28 . As 
expected, AIMD-Z overestimates T m , and the amount of the overestimate T\ lq — T m is 
similar to our prediction from eqn 
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4 Discussion and conclusions 



Our investigation has yielded several simple but important insights into the kinetics of 
homogeneous melting under (N,V,E) conditions. First, we have confirmed the existence 
of a "limit of superheating" X"ls proposed in previous work [7J, and we have shown that 
for a given "excess" quasi-steady temperature T ao \ — Tls of the initial superheated solid, 
and for a given number of atoms A, there is a fairly well defined probability per unit 
time (reciprocal of mean waiting time (r w )) of making the melting transition. For given 
N, (t w ) appears to scale roughly as A/{T\\ q — T m ) 2 with the excess of the final liquid 
temperature Ti; q above the true thermodynamic melting temperature T m associated with 
the specified liquid density. The coefficient A appears to be roughly proportional to 1/A. 
These rather simple findings are clearly interesting, and it is natural to ask whether they 
will hold true for materials other than the particular transition metal studied here. At 
present, we have no way of answering this question, and there is now a clear need to 
extend the investigation to materials of other kinds. 

Another natural question concerns the relation between homogeneous melting and 
the kinds of metastablc behaviour well known in, for example, supersaturated vapours or 
supercooled liquids. Of course, theories of such phenomena have an extremely long history, 
the point of reference often being "classical nucleation theory" (CNT) [45 . In CNT, the 
waiting time for the transition to the thermal equilibrium state (condensation, freezing,...) 
is governed by the time needed to form a "critical nucleus"; the associated free-energy 
barrier results from a competition between the lowering of bulk free energy resulting from 
the transition and the free energy increase due to the formation of interfaces (liquid- 
vapour, solid-liquid,...) during the transition. The metastable behaviour of a superheated 
solid that is only slightly beyond the thermal-equilibrium stability field of the solid should 
also be describable by an appropriate CNT. However, it is possible that the formation of 
a critical nucleus described by such a CNT has little to do with the homogeneous melting 
observed in previous simulations [B] and in the ones presented here, for two reasons. 
First, CNT theories and other theories of nucleation would not predict a "superheating 
limit" . Instead, they would predict a mean waiting time for nucleation that decreases 
continuously as we move further into the field of thermodynamic instability. Second, 
the transition associated with classical nucleation will generally lead to a final state in 
which both phases are present, rather than the single (liquid) phase seen in simulations 
of homogeneous melting. 

The persistent alternation between solid and liquid states that we observe for the 
very small 96-atom system is relevant here. To understand why this happens, and to 
know when we should expect to see it, we recall the ergodic principle, which is generally 
accepted to hold for condensed-matter systems. This states that the trajectory produced 
by (A, V, E) m.d. starting from any phase-space point (set of positions and momenta) 
having specified total energy E will, if continued long enough, pass arbitrarily close to an 
arbitrarily chosen phase-space point of the same E. The configurations we are concerned 
with here are either solid-like or liquid-like: in none of the simulations we have performed 
do we see solid and liqid simultaneously present, so that the configurations are either 
one or the other, except for the very small fraction of configurations that occur during 
the transitions from solid to liquid or vice versa. Now suppose that, starting from the 
superheated solid, the system has homogeneously melted and become liquid. Then the 
ergodic principle tells us that the given trajectory, if continued long enough, will eventually 
re-enter regions of solid-like configurations, so that it will re-freeze. Indeed, the trajectory 
will densely cover the entire constant-i? manifold, and will spend well defined fractions 
aiiq and a so i = 1 — an q in the liquid and solid states. This is exactly what we have seen, 
and Fig. [7] displays the value of au q as a function of 7ii q . 

It is straightforward to confirm that this explanation is correct. The fractions an q and 
a so i are proportional to the numbers Wy lq and W so \ of liquid-like and solid-like microstates 
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on the given constant-T manifold. But these are related to the entropies Sn q and S so \ of 
the corresponding macrostates: S B ol,Hq = &b In Ws i,ii q . Hence, we have: 

"uq/asoi = aiiq/(l - aiiq) = cxp((5n q - 5 so i) /fc B ) • (3) 

This means that if we choose E = E so that the system spends equal amounts of its 
time in solid and liquid states, then the entropies are equal. Let the temperatures in this 
situation be T s ° ol and Xjj q . If we go to a nearby energy E = E + SE, then an q and a so i 
change, because: 



But (dS/dE)v = 1/T, and SE can be expressed as SE ~ C v ,iiq5T nq , where STu q = 
T\i q — Tj? q , and C v ,ii q is the isochoric specific heat of the liquid. From eqns |3| and Q, 
we then have: 

aiiq = i + cx P (-<jT liq /r int ) ' (5) 

where the temperature interval Tj nt over which the transition occurs is: 

_L = (iVC v ,Uq/feB)(^-i) ■ (6) 

Jmt \ liq sol / 

The Fermi function of eqn ^ is the function we have used to fit our simulation results 
for oey lq as a function of T nq (Fig. [7]), and the parameters that emerge from this fit are 
= 6705 K, Ti n t = 124 K. We can now check that this value of T- mt obtained by empirical 
fitting is indeed consistent with the prediction of eqn From our EAM simulations of 
the liquid, we obtain the estimate C V) n q /fcB = 3.36. Using the observed values of T® ol and 
T^, we then obtain the prediction T nt = 125 K, which is very close (perhaps fortuitously 
close) to what we obtain from fitting. 

It is clear from what we have said that solid-liquid alternation is only an important 
effect for small systems. The key feature that makes it easy to observe in our 96-atom 
system is that the temperature distributions of the solid and liquid states overlap signif- 
icantly (Fig. [6| . Since the rms temperature fluctuation of a single-phase system in the 
microcanonical ensemble is proportional to 1/yN, the overlap becomes negligible for large 
systems. The same conclusion is clear from eqns ^ and ([6]), which show that temperature 
interval T nt goes as 1/N. For a large system, once the quasi-steady temperature T so \ of 
the initial solid exceeds Tls , the homogeneous melting transition is effectively irreversible. 

Our results shed light on the Z method for the determination of melting properties. 
This method is simple to use, but our work shows that it generally gives only an upper 
bound to the melting temperature associated with a given liquid density, unless measures 
are taken to correct it. This is because melting may not be observed even when the quasi- 
steady temperature T so \ of the solid is above Tls- Indeed, melting is very unlikely to be 
seen if T so \ — Tls is such that (r w ) is much longer than the duration of the simulation. 
This is a particular problem for AIMD, where we have shown for the case of hep Fe that, 
even with 50 ps simulations on a system of 150 atoms, melting is unlikely to be seen until 
T so i — Tls — 300 K, in which case the final liquid temperature will overestimate T m by 
~ 300 K. For the systems of less than 100 atoms and simulations of less than 10 ps used 
in some recent AIMD-Z work [32j[33], the overestimation is likely be much worse. 

It is clearly desirable to have ways of correcting for the overestimate of T m given by 
the Z method. Our work demonstrates that calculation of the mean waiting time (r w ) 
provides one way of doing this. For large enough T so \ — Tls, melting will occur rapidly, and 
(r w ) can then be estimated by repeating the simulation many times so as to reduce the 
statistical errors on (t w ). If this is done at two or more values of T so \ — Tls, we then have 
information about the dependence of (t w ) -1 / 2 on Ti q , from which the necessary correction 
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can be made. This would be a somewhat expensive procedure for AIMD, but would have 
the great advantage of being simple and automatic, since many simulations could be run 
simultaneously on a large compute cluster. The relation with parallel replica methods 46 
will be noted. As an alternative, it may well be possible to use Bayesian techniques to 
extract the information needed for corrections from the sequences of simulations that 
are required in any case by the Z method. On a completely different line of thought, 
we remark that since homogeneous melting is a rare-event problem, it may be possible 
to exploit techniques used for other rare-event problems to accelerate melting in the 
Z method. Metadynamics [UJ might be one such technique, since it would be easy to 
adopt "collective variables" that would discourage the system from remaining too long in 
the solid state [H]. We plan to investigate some of these possibilities in the future. 

In conclusion, our m.d. simulations on the homogeneous melting of the transition 
metal Fe confirm the existence of a rather well defined limit of superheating, beyond 
which melting occurs on a typical time-scale of ns or less. We have shown that the 
statistical distribution of waiting times r w before melting displays a typical "rare event" 
character, consistent with a probability per unit time for melting to occur. The mean 
waiting time (t w ) lengthens rapidly as the superheating limit is approached from above, 
being roughly proportional to the inverse square of the excess beyond the limit; it also 
lengthens as the system size (number of atoms N) decreases, being roughly proportional to 
1/N. This means that the Z method for calculating melting temperatures can be subject 
to large errors if it is applied to small systems over short times, though the method can 
work successfully under suitable conditions. We have noted that these conditions have 
not always been satisfied in earlier work. 
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Figure 1: Time-dependent temperature and pressure in four independent simulation runs, 
showing homogeneous melting from superheated hep solid Fe in a system of 7776 atoms. All 
four simulations were initiated from perfect crystal positions, with initial random velocities 
corresponding to the same temperature T m = 15600 K, the mean quasi-steady temperatures 
of the superheated solid and the final liquid being T so \ = 7590 K and Tu q = 6315 K. 
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Figure 2: Histograms of waiting times r w before the transition to liquid constructed from 
repeated simulations at two initial temperatures for the system of 7776 atoms. Histograms 
shown by dashed (red) and solid (black) lines result from initial temperatures of 2] = 15800 
and 16000 K respectively, the quasi-steady solid and liquid temperatures in the two cases 
being T so \ = 7640 and 7740 K and 7i iq = 6410 and 6505 K. Dashed and dotted curves show 
exponential functions fitted to histograms (see text). 
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Figure 3: Dependence of mean waiting time (t w ) on final liquid temperature T\i q for systems 
of N = 7776 (black circles), 976 (red squares), 392 (green diamonds), 150 (blue triangles) 
and 96 (brown stars) atoms. Quantity plotted is (tw)^ 1 / 2 as function of Tu q . Straight lines 
are linear least-squares fits to the data for each N value. 
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Figure 5: Alternation between solid and liquid states: temperature as function of time in 
one of the constant-energy m.d. simulations on the system of 96 atoms, with total energy 
such that mean liquid-state temperature is Tn q = 6760 K, showing alternation between mean 
temperatures T so \ and Tj iq . 
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Figure 6: Histograms of temperature distribution at different constant total energies E in the 
system of 96 atoms. The histogram at each E was obtained by sampling over typically 128 
simulations, each having a typical duration of 5 ns. Instead of giving E directly, we specify 
each histogram by the liquid-state temperature Ti; q . Histograms shown by solid (black), 
dashed (red), dotted (green), chain (blue) and dotted-chain (black) lines are for Tu q = 6935, 
6760, 6590, 6565 and 6473 K. 
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Figure 7: Fraction of time spent by the system in the liquid state for different liquid-state 
temperatures 71i q in simulations of 96-atom system. 
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Figure 8: Z plot from a sequence of constant-energy AIMD simulations on system of 150 
atoms, duration of simulations being 50 ps. Black filled circles with error bars show final 
mean temperature and pressure, with upper-left branch corresponding to energies for which 
the system remains solid, and lower right branch to energies for which homogeneous melting 
occurs. Dashed (red) line shows the ab initio melting curve obtained in earlier work using 
the free energy technique, and green filled square with error bar shows point on ab initio 
melting curve obtained with ab initio m.d. simulations on a system of 980 atoms containing 
coexisting solid aid liquid. 
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